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Abstract. The analytical expression for integral energy spec- 
tra and zenith angle distributions of atmospheric muons at 
large depths is derived. Fluctuations of muon energy losses 
are described using the parametrized correction factor. The 
fitting formula for the sea level muon spectrum at different 
zenith angles for spherical atmosphere is proposed. The con- 
crete calculations for pure water are presented. 



1 Introduction 

The last several years have been marked by the start of 
full scale data taking of large neutrino and muon tele- 
scopes located at lake Baikal (NT-200) and in deep polar 
ice (AMANDA). Two underwater telescopes ANTARES and 
NESTOR assuming installation at greater depths are being 
intensively constructed in the Mediterranean. The possibiUty 
of deployment of telescopes with huge detecting volumes up 
to 1 km'^ is also under wide investigation. 

So, the knowledge of expected angular distribution of inte- 
gral flux of atmospheric muons deep underwater is of interest 
not only for cosmic ray physics but also for the estimation of 
the possible background for neutrino detection and at last for 
a test of the correctness of underwater telescope data inter- 
pretation by using the natural flux of atmospheric muons as 
calibration source. The last item frequently implies the es- 
timation with an appropriate accuracy (e.g., better than 5 % 
for a given sea level spectrum) the underwater integral muon 
flux for various sets of depths, cut-off energies and angular 
bins especially for telescopes of big spacial dimensions. 

Up to now the presentation of the results of calculations 
of muon propagation through thick layers of water both for 
parent muon sea level spectra (especially for angular depen- 
dence taking into account the sphericity of atmosphere) and 
for underwater angular flux has not been quite convenient 
when applied to concrete underwater arrays. In addition, a 
part of numerical results is available only in data tables (of- 



ten insufficient for accurate interpolation) and figures. The 
possibility of direct implementation of Monte Carlo methods 
depends on the availability of corresponding codes and usu- 
ally assumes rather long computations and accurate choice 
of the grid for simulation parameters to avoid big systematic 
errors. Therefore, there remains the necessity of analytical 
expressions for underwater muon integral flux. In addition, 
the possibility of reconstructing the parameters of a sea level 
spectrum by fitting measured underwater flux in the case of 
their direct relation looks rather attractive. 

In this paper we present rather simple method allowing one 
to analytically calculate the angular distribution of integral 
muon flux deep under water for cut off energies (1-10^) GeV 
and slant depths of (1-16) km for conventional (tt, K) sea 
level atmospheric muon spectra fitted by means of five pa- 
rameters. The fluctuations of muon energy losses are taken 
into account. 



2 Basic formulas 

According to the approach of work (Klimushin, Bugaev, and 
Sokalski, 2000) the analytical expression for calculations of 
underwater angular integral flux above cut-off energy Ef for 
a slant depth R — h/ cos 9 seen at vertical depth h at zenith 
angle 6 and allowing for the fluctuations of energy loss is 
based on the relation 



Ffi{>Ef,R,i 



Fci{>Ef,R,e) 
Cf{>Ef,R,0)' 



(1) 
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where correction factor C/ is expressed, by definition, by 
the ratio of theoretical integral flux calculated in the contin- 
uous loss approximation to that calculated by exact Monte 
Carlo. The analytical parametrization for C/ is presented in 
Refs. (Klimushin, Bugaev, and Sokalski, 2000, 2001). The 
dependencies of the correction factor on Ef and R, calcu- 
lated for any reasonable sea level spectrum represent the set 
of rather smooth curves. 

The angular flux Fci{> Ef, R, 9) based on effective linear 
continuous energy losses a + (3E having 2 slopes, is calcu- 
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lated by the following rule: 



F^l{>Ef,R,e;au0i) 
Fci{>Ei2,iR-Ri2),e; a2,/?2) 



Fei{>Ef,R,0) = 
fori? < Ri2, 



fori? > Ri2- 



(2) 



Here E12 is the energy in the point of slope change from 
(ai, /3i) to {a2,P2) and R12 is the muon path from the en- 
ergy £^12 till Ef which is given by 



i?i2 = In 
Pi 



ai + E12P1 
ai + Ef0i 



The formula for integral muon angular flux in the assump- 
tion of Unear continuous energy losses is as follows: 



F^i{>Ef,R,0; a,f3) 



7 



X ^ Do,Ei:{e){Es+yi)-^{l-Zif-'^S{zi,^), (3) 

i—-K,K 

where subscript i stands over both pion (tt) and kaon {K) 
terms and 

_ Es:{9)e-^^ _ Es:m 

Zi - — — : , Eq. [U) — 



Ef+Vi ' 



cos6* 



00 ^ n s — 1 

s{zn) = i + Y.^\z-[\{{.i + 3)) =i + rTT 



2z' 



6z-^ 



(7 + l)(7 + 2) (7 + l)(7 + 2)(7 + 3) 

When using expression (3) for slant depths R > R12 one 
must substitute R ^ (R — R12 ) and Ef E12 and use the 
values (a2 , P2) for a loss description. For slant depths R < 
R12 the use of (3) remains unchangeable and the loss values 
are expressed by (ai , /3i). This algorithm may be extended to 
computations with any number of slopes of the energy losses. 



those of the differential sea level muon spectrum, for which 
we use the following parametrization: 



D{Eo,9) 



^0 



E ■ 

i=n,K 



Oi 



+ Eo/E^lie) ' 



(4) 



where 7 is a spectral index and Eq^ ^ (9) have approximate 
sense of critical energies of pions and kaons for given zenith 
angle and Eq"^ ^ (0°) are those for vertical direction. The cor- 
responding angular distrubution should be introduced using 
an analytical description of effective cosine cos 9* taking into 
account the sphericity of atmosphere. It should be noted that 
the description of underwater angular flux with the 5 parame- 
ters of a sea level spectrum gives the possibility of their direct 
best fit by using the experimental underwater distribution. 

Flux value in (3) is expressed in units of (cm~^s~^sr~^) 
and all energies are in (GeV), slant depth R in units of 
(gcm~^), loss terms a and (3 in units of (10~"^GeVcm^g~^) 




E„ (GeV) 



Fig. 1. Differential spectra of conventional muons at sea level for 
six zenith angles, sec 6: 1, 2, 3, 5, 10, and 20, from bottom to top. 
Curve labels correspond to values of sec d. Solid curves: spectrum 

based on data tables from (Misaki et al, 1999), dotted curves: spec- 
trum defined by Eq. (5), dashed curves: spectrum of Ref. (Volkova, 
Zatsepin, and Kuz'michev, 1979) shown only for two values of 
sec 9: 1 and 20. 

and (10~^cm^g^^), correspondingly. For the description of 
effective linear continuous energy losses we use the follow- 
ing values of parameters when substituting in (2): (q:i=2.67, 
^1=3.40) and (a2=-6.5, /32=3.66) with Si2=35.3 TeV. 

To examine the angular behaviour of a flux given by the 
formula (1) by means of the comparison with numerical cal- 
culations we used the following parameters of the sea level 
muon spectrum: 

-3 



The 5 parameters (i?o, , ^o^ , ^SI (0°), (0°), 7) are (0°) 



£>o. = 0.175, 
103 GeV, E^l{0°) 



Doj, = 6.475 X 10" 
810 GeV, 7 = 2.72 . 



These values have been chosen according to splines com- 
puted in this work via the data tables kindly given us by au- 
thors of Ref. (Misaki a/., 1999). When checking the values 
of fit spectrum for cos 6'=(0.05-l .0) we realized that the stan- 
dard description of effective cosine (with geometry of spher- 
ical atmosphere and with definite value of effective height of 
muon generation) is not enough and one should introduce an 
additional correction S{9) leading to (10-20) % increase of 
effective cosine value for cos^ < 0.1. The reason of an ap- 
pearing of this correction is that the concept of an effective 
generation height is approximate one. It fails at large zenith 
angles where the real geometrical size of the generation re- 
gion becomes very large. 

We should note that our expression (4) for the sea level 
muon spectrum does not contain a contrubution from atmo- 
spheric prompt muons. According to the most recent calcula- 
tions based on perturbative QCD, this contribution becomes 
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essential only at Eq > 10® GeV. The predictions of non- 
perturbative models are slightly more optimistic. We plan 
to generaUze our approach and include this contribution in 
our following paper. Incidentally, inclusion of prompt muons 
should be done in parallel with taking into account the steep- 
ening of the sea level muon spectrum due to the knee in the 
primary cosmic ray spectrum. 

The resulting fit of angular sea level spectrum in units 
of (cm~^s~^sr~^Get^~^) is given by 



1 



Eocos6»** 



+ 



0.037 



1 + 



i;ocos6'** 



1 + 



103 810 
with modified effective cosine expressed by 

cosr* = 5(61) cos r, 



(5) 



(6) 



where cos 9* is derived from spherical atmosphere geometry 
and is given by the polynomial fit: 



cos 61* 



i=0 



i cos* 9, 



(7) 



with the coefficients of the decomposition assembled in Ta- 
ble 1. The accuracy of (7) is much better than 0.3 % except 
the region cos ^=(0.3-0.38) where it may reach the value of 
0.7 %. Note that for cos 9 > 0.4 the influence of the curva- 
ture of real atmosphere is less than 4 % but for cos ^ < 0.1 it 
is greater than 40 %. 

S{9) is the correction which is given for sec ^ < 20 by 



S{9) = 0.986 + 0.014 800 6*. 



(8) 



Fig. 1 illustrates the limits of applicability of angular spec- 
trum given by Eq. (5), for energy and zenith angle variables. 
The energy region, inside which the deviation from parent 
spectrum is less than 5 %, is shifted from (0.3-200) TeV for 
cos 61=1.0 to (1.5-300) TeV for cos 61=0.05. The sea level 
spectrum given by (5) is valid only below the knee {Eq ^ 300 
TeV) of primary cosmic ray spectrum. 

Correspondingly, for critical energies in expression (3) one 
should use cos 9** instead of cos 9*. 



3 Comparison with numerical calculations 

The examination of (3) showed rather quick convergence of 
series 8(2,7) with increase of R and Ef. Therefore, for 
the accuracy of computation better than 0.1 % it is quite 
enough to take only four first terms of this series (up to z^) 
for all values R > I km and Ef in (1-10'') GeV. Even us- 
ing the two terms leads to the accuracy of 1.3 % for (i?=1.15 
km, Ef=l GeV) and <0.5 % for (R > 2.5 km, Ef > 1 GeV). 
Fig. 2 shows the comparison of underwater angular integral 
fluxes aUowing for loss fluctuations at different basic depths 
h (of location of existing and planned telescopes) calculated 
both numerically using MUM code (Sokalski, Bugaev, and 
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Fig. 2. Underwater integral muon flux aUowing for loss fluctuations 
as a function of zenith angle at different vertical depths, for cut- 
off energy Ef =10 GeV. Four curves correspond to vertical depths 
h: 1.15 km, 1.61 km, 2.0 km, and 3.0 km, from top to bottom. 
Solid curves result from numerical computations by using the sea 
level spectrum based on data tables from (Misaki et al, 1999) and 
MUM code of muon propagation. The dots are result of analytical 
calculation with expression (1) by using the sea level spectrum (5). 



Klimushin, 2000) for parent sea level spectrum and analyti- 
cally (1) for the spectrum given by (5). 

We realized that the error given by formula (1) for all men- 
tioned sea level spectra is within the corridor of ±(4-6) % for 
all cut-off energies E'/=(l-10^) GeV and slant depths i?=(l- 
16) km (corresponding angle is expressed by 0086* ~ h/R 
for a given vertical depth h). This is proved for h 'ma range 
(1-3) km. For bigger cut-offs of i;/=(l-10)TeV the corri- 
dor of errors is ±(5-7) % for i?,=(l-13) km. Note that for the 
sea level spectrum (5), just used for C/ parametrization, the 
errors are smaller on 2 %. 

The accuracy of the parametrization, used for the correc- 
tion factor as a function of Ef and slant depth R is rather 
high and is about ±5 % for all angles and kinds of the sea 
level spectrum (assuming that the spectral index 7 is approx- 
imately within (2.65-2.78)). It results in the possibiUty to 
use it for an estimating numerically from various sea level 
spectra the value of an angular integral flux allowing for fluc- 
tuations of losses without direct Monte Carlo simulations. 

Note that the expression (1) may be directly used for an 
ice after substitution R^ R/ p, with p being the ice density, 
and, with an additional error of ~ 2 %, for a sea water. In 
spite of seeming complexity of the formulae (1), (2) and (3) 
they may be easily programmed. 

The validity of proposed formula up to cut-off energies 
10 TeV allows a calculation of underwater angular differen- 
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Table 1. Coefficients d of the fitting formula (7) for effective cosine witli the maximum relative errors. 
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Fig. 3. Underwater differential muon spectrum allowing for loss 
fluctuations as a function of energy Ef at vertical depth of 1.15 km. 

The curves are given for twelve zenith angles cos6: 1.0, 0.9, 0.8, 
0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.15, 0.1, and 0.05, from top to bot- 
tom. Solid curves result from numerical differentiation of integral 
flux computed by using the sea level spectrum based on data tables 
from (Misaki et al., 1999) and MUM code of muon propagation. 
Dotted curves result from numerical differentiation of the analyti- 
cal expression (1) based on the sea level spectrum (5). 

tial spectrum D{Ef,R,9) by means of numerical differenti- 
ation of expression (1). It leads to rather appropriate results 
up to slant depths (11-12) km. We illustrate this in Fig. 3 
by comparison the underwater angular spectra calculated by 
numerical differentiation of integral flux computed by using 
the sea level spectrum based on data tables from (Misaki et 
al., 1999) and MUM code of muon propagation, and numer- 
ical differentiation of the analytical expression (1) based on 
the sea level spectrum (5). For a vertical depth /i=1.15 km 
it leads to errors ±4% for Ef =(20-8 x 103)GeV for tiie 
angles corresponding to slant depths R ~ h/ cos 9 of (1- 
3) km and ±(6-8) % for Ef =(30-5 x lO^) GeV for flie slant 
depths (3-12) km. Even for i?=23.2 km the result is stiU valid 
within ± 1 % but for the very narrow energy region E / =(90- 
300) GeV. 



4 Conclusions 

The analytical expression presented in this work allows to es- 
timate for fluctuating losses the integral flux of atmospheric 
muons in pure water expected for different zenith angles, 
cos 6'=(0. 05-1.0), at various vertical depths at least of h=(l- 
3) km for different parametrizations of the sea level muon 
spectra. The errors of this expression are estimated to be 
smaller than ±(4-6) % for cut-off energies ranged in Ef=(l- 
10^) GeV and slant depths inh/ cos 9=(l-l6) km. The main 
advantage of the presented formula consists in the possibility 
of the direct best fit of 5 parameters of parent sea level spec- 
trum using angular distribution of underwater integral flux 
measured experimentally at a given vertical depth. 

The validity of this analytical expression with an accuracy 
of ±(5-7) % for Ef=(lO^-lO'^) GeV and slant depths of (1- 
12) km gives also the possibility of estimation the angular 
underwater differential spectrum (by means of numerical dif- 
ferentiation) with error smaller than ±(6-8) % for energies of 
(30-5 xlO^) GeV. 

The accuracy of the parametrization, used for the correc- 
tion factor as a function of and slant depth R is rather 
high and is about ±5 % for all angles and kinds of the sea 
level spectrum (assuming that the spectral index 7 is approx- 
imately within (2.65-2.78)). It results in the possibihty to 
use it for an estimating numerically from various sea level 
spectra the value of an angular integral flux allowing for fluc- 
tuations of losses without direct Monte Carlo simulations. 

The proposed method may be adapted to estimations in 
rock after corresponding description of the correction factor 
and continuous effective losses. 
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